     program read_pyscf
!
     use mod_ri
     use mod_ioff
     use mod_iop
     use mod_orbit
     use mod_size
     use mod_symm
     implicit none
!
     integer irrepi, i, j, k, irrepna, irrepnb
     integer isymoff(8,8,25), irrepa2(255,2), irpdpd(8,22)
     integer IOFF(2), IFLAGS(100), nsize
     integer n, nmo, symmo(8), istat, memsize
     integer i0, i1, i2, i3, i4, i5, i6, i7, i8, i9
     integer i10, i11, i12, i13, i14, i15, i16, i17, i18, i19
     integer irrepci, ioffr0, ioffs0, ioffr, ioffs, numi
     integer idx, irrep, irrept, nvmax, nomax, idxn
     integer l1, l2, mgb, sym_id(8), frozen
     integer irrepcd, irrepab, numab, irrepd, irrepc, numc, numd
     integer ioffin, ioffout1, ioffout2, c, d, iin, iout1, iout2
     integer ioffabc(9), ioffs2, ioffs8(8), numci
     integer irrepa, numa, a, ain
     real*8 escf, rt1, rt2

     integer, allocatable :: isymmo(:,:), ioffmo(:,:), irrepmo(:,:)
     integer, allocatable :: isymmo_u(:,:,:)
     integer, allocatable :: isymvrta(:,:), isymvrtb(:,:)
     integer, allocatable :: isympopa(:,:), isympopb(:,:)
     integer, external :: idsymsoc, irpdso, isymoffso
     real*8, allocatable :: ea(:), eb(:), ea_resort(:), eb_resort(:)
     real*8, allocatable :: faiaa(:), faibb(:) 
     real*8, external :: nrm2so
     real*4, allocatable :: icore(:)
     real*4, allocatable :: itmps(:)
     character*3 :: group, sym_name(8)
     character*12 :: inametmp
     character*12 :: string2
     character*30 :: iname(100)
     logical yesno
!---fr mod_iop
!
     write(6,*) 'Call cc_pyscf'
     prec=8; isd=2; igpu=0; totgpu=0
     call readinp
     igpu=0; totgpu=0
!
     write(6,*)'finished readinp'
     isd=prec/4    
     frozen = 0 
!
     cc_int = 2  ! for choose where to get the integrals, =1 for cfour and =2 for pyscf
!
     open(unit=1100,file="cc_pyscf_int")
     if(ldebug==.true.) then 
        write(6,*)'finished open cc_pyscf_int'
     endif
     pop=0; vrt=0
     read(1100,*) noa, nob, nva, nvb, nmo, escf, group
!    pop(1,1)=noa; pop(1,2)=nob
!    vrt(1,1)=nva; vrt(1,2)=nvb
     if(ldebug==.true.) then 
         write(6,*)'nv and no', nva, nvb, noa, nob, nmo, escf
     endif
     nirrep=1
     allocate(ea(noa+nva), eb(nob+nvb), stat=istat)
     allocate(ea_resort(nmo), eb_resort(nmo), stat=istat)
     allocate(itmps(noa*noa*nva))
!     call writecommon
     read(1100,*)(ea(i), i=1, noa+nva)
     read(1100,*)(eb(i), i=1, nob+nvb)
!-----------------------------------------------------------
!-----------------------------------------------------------
     isymv = 1

     allocate(isymmo(nmo, 8), isymmo_u(nmo,8,2))
     write(6,*)'iop_ref', iop_ref
!--------read symm and orbit in "mod_symm and mod_orbit"
     if(iop_ref==0) then 
        open(22, file='SYMM', form='formatted')
        rewind 22
        read(22,*) nmo
        read(22,*) nirrep
        read(22,*) dirprd
        read(22,*) vrt(:,1)
        read(22,*) pop(:,1)
        read(22,*) symmo
        read(22,*) isymmo
        read(22,*) symrx, symry, symrz 
        close(22)
     elseif(iop_ref==1) then 
        open(22, file='SYMM_UHF', form='formatted')
        rewind 22
        read(22,*) nmo
        read(22,*) nirrep
        read(22,*) dirprd
        read(22,*) vrt
        read(22,*) pop
        read(22,*) symmo
        read(22,*) isymmo_u
        read(22,*) symrx, symry, symrz
        close(22)
        isymmo = isymmo_u(:,:,1)
     endif
     

     do i = 1, nirrep
        read(1100,*) sym_name(i), sym_id(i)
     enddo
     close(1100)
     if(iop_ref==0) then
        do i = 1, nirrep  
           vrt(i,2) = vrt(i,1) 
           pop(i,2) = pop(i,1)
        enddo
     endif
!     write(6,*)'vrt', vrt
     vrta = vrt(:,1); vrtb = vrt(:,2)
     popa = pop(:,1); popb = pop(:,2)

     allocate(ioffmo(nmo,2), irrepmo(nmo,2))
     idx = 0
     do n = 1, nmo
        do irrep = 1, nirrep
        do i = 1, pop(irrep,1)+vrt(irrep,1)
           if(isymmo(i, irrep)==n) then
              idx = idx + 1
              irrepmo(n,1) = irrep
              ioffmo(n,1) = i
           endif
        enddo
        enddo
     enddo

     if(iop_ref==1) then 
       idx = 0
       do n = 1, nmo
          do irrep = 1, nirrep
          do i = 1, pop(irrep,1)+vrt(irrep,1)
             if(isymmo_u(i, irrep,2)==n) then
                idx = idx + 1
                irrepmo(n,2) = irrep
                ioffmo(n,2) = i
             endif
          enddo
          enddo
       enddo
     endif  
     write(6,*)'sym_name', sym_name, noa+nva
!     frozen = 11 
     write(6,*)'==========alpha orbital energy================='
     do i = 1, noa
        irrep = irrepmo(i,1)
!        write(6,*)'irrep', irrep
        write(6,'(T10, I6, T20, F18.8, T45, I1, T50, A12)')&
                    i+frozen, ea(i), irrep, sym_name(irrep)
     enddo
     write(6,*)'==============================================='
     do i = noa+1, noa+nva
        irrep = irrepmo(i,1)
!        write(6,*)'irrep', irrep
        write(6,'(T10, I6, T20, F18.8, T45, I1, T50, A12)')&
                    i+frozen, ea(i), irrep, sym_name(irrep)
     enddo
     if(iop_ref==1) then 
     write(6,*)'===========beta orbital energy================='
     do i = 1, nob
        irrep = irrepmo(i,2)
!        write(6,*)'irrep', irrep
        write(6,'(T10, I6, T20, F18.8, T45, I1, T50, A12)')&
                    i+frozen, eb(i), irrep, sym_name(irrep)
     enddo
     write(6,*)'==============================================='
     do i = nob+1, nob+nvb
        irrep = irrepmo(i,2)
!        write(6,*)'irrep', irrep
        write(6,'(T10, I6, T20, F18.8, T45, I1, T50, A12)')&
                    i+frozen, eb(i), irrep, sym_name(irrep)
     enddo
     endif

     idx = 0 
     do irrep = 1, nirrep
        do i = 1, popa(irrep) 
           idx = idx+1
           ea_resort(idx) = ea(isymmo(i, irrep))
        enddo
     enddo 

     do irrep = 1, nirrep
        do i = popa(irrep)+1, popa(irrep)+vrta(irrep)
           idx = idx+1
           ea_resort(idx) = ea(isymmo(i, irrep))
        enddo
     enddo

   if(iop_ref==1) then 
     idx = 0
     do irrep = 1, nirrep
        do i = 1, popb(irrep)
           idx = idx+1
           eb_resort(idx) = eb(isymmo_u(i, irrep,2))
        enddo
     enddo

     do irrep = 1, nirrep
        do i = popb(irrep)+1, popb(irrep)+vrtb(irrep)
           idx = idx+1
           eb_resort(idx) = eb(isymmo_u(i, irrep,2))
        enddo
     enddo
    elseif(iop_ref==0) then
       eb_resort = ea_resort
    endif 

     do i = 1, noa+nva
!         write(6,'(A, I8, 2F18.12)')'ea and eb', i, ea_resort(i), ea(i)
     enddo
!
!====save energy for CC calculation
     if(prec==8) then       
        if(iop_ref==0) then 
          open(unit=77, file='energydatar', form='unformatted')
          rewind 77
          write(77)(EA_resort(i), i = 1, noa+nva)
          write(77) ESCF
          close(77)
        elseif(iop_ref==1) then 
          open(unit=77, file='energydatau', form='unformatted')
          rewind 77
          write(77)(EA_resort(i), i = 1, noa+nva)
          write(77)(EB_resort(i), i = 1, noa+nva)
          write(77) ESCF
          close(77) 
        endif 
     elseif(prec==4) then 
        if(iop_ref==0) then 
           open(unit=77, file='energydatars', form='unformatted')
           rewind 77
           call dtos(EA_resort, itmps, noa+nva)
           write(77)(itmps(i), i = 1, noa+nva)
           write(77) ESCF
           close(77)
        elseif(iop_ref==1) then 
           open(unit=77, file='energydataus', form='unformatted')
           rewind 77
           call dtos(EA_resort, itmps, noa+nva)
           write(77)(itmps(i), i = 1, noa+nva)
           call dtos(EB_resort, itmps, noa+nva)
           write(77)(itmps(i), i = 1, noa+nva)
           write(77) ESCF
           close(77) 
        endif 
     endif
     nvmax = 1; nomax = 1
     do i = 1, 8
        do j = 1, 2 
           nvmax = max(vrt(i,j), nvmax)
           nomax = max(pop(i,j), nomax)
        enddo
     enddo
!    do i = 1, 8
!       write(6,'(A, 8I4)')'dirprd', dirprd(:,i)
!    enddo

     nva = sum(vrta)
     noa = sum(popa)
     zorb(1)=1
     zorb(2:8)=0

!  sae other information for cfour calculation 
     call writecommon
!-----------------------------------------------------------
!-----------------------------------------------------------
!---fr mod_size and mod_ioff
!    mgb=40
!    icrsiz = 1024*1024*1024/8*40  !80GB   memory size 
!    igpusize= 1024*1024*1024/8*10 !10GB   gpu size
!    CALL READINP
!     icrsiz = 1024*1024*1024/8*mgb  !80GB   memory size 
 
     allocate(icore(icrsiz), STAT = istat)      
     call checkalloc(istat, 'icore')
!
     if(ldebug==.true.) then 
        write(6,*)'iprcln', iprcln, iprcwd, ichcsz
        write(6,*)'sym', nirrep
        write(6,'(A,8I4)')'vrta', vrta 
        write(6,'(A,8I4)')'popa', popa
        write(6,'(A,8I4)')'vrtb', vrtb 
        write(6,'(A,8I4)')'popb', popb
        do i = 1, 8 
           write(6,'(A,8I4)')'DIRPRD', DIRPRD(:,i)
        enddo
!        write(6,*)'vrt', vrt
     endif
     call walltimeso(rt1)
     nvmax = 1; nomax = 1
     do i = 1, 8
        do j = 1, 2 
           nvmax = max(vrt(i,j), nvmax)
           nomax = max(pop(i,j), nomax)
        enddo
     enddo
     allocate(isymvrta(nvmax, 8), isymvrtb(nvmax,8), isympopa(nomax,8),&
                         isympopb(nomax,8), STAT = istat)
     call checkalloc(istat, 'isymvrtra ...')
     isymvrta = 0; isymvrtb = 0; isympopa =0; isympopb = 0
     do irrept = 1, nirrep
        do i = 1, vrta(irrept) 
           isymvrta(i,irrept) = isymmo(popa(irrept)+i, irrept)
        enddo
        do i = 1, popa(irrept) 
           isympopa(i,irrept) = isymmo(i, irrept)
        enddo
        if(iop_ref==0) then 
           isymvrtb = isymvrta
           isympopb = isympopa
        elseif(iop_ref==1) then 
           do i = 1, vrtb(irrept) 
              isymvrtb(i,irrept) = isymmo_u(popb(irrept)+i, irrept,2)
           enddo
           do i = 1, popb(irrept) 
              isympopb(i,irrept) = isymmo_u(i, irrept,2)
           enddo
        endif
     enddo  
     if(ldebug==.true.) then 
        write(6,*)'precision in cc_pyscf_rhf', prec
     endif
!=====resort fae 
!=====
     call fai_resort(icore, isympopa, isympopb,  &
                isymvrta, isymvrtb, nvmax, nomax)
!======================
     call flush(6)
     call walltimeso(rt1)
     if(iop_ref==0) then 
         write(6,*)'close shell'
         if(prec==8.and.l_ri==.false.) then 
         call run_rhf(icore, isympopa, isympopb, isymvrta, isymvrtb, isymmo, &
                                      nvmax, nomax, nmo) 
         elseif(prec==8.and.l_ri==.true.) then 
         call run_rhf_ri(icore, isympopa, isympopb, isymvrta, isymvrtb, isymmo, &
                                      nvmax, nomax, nmo, icrsiz/isd)
         elseif(prec==4.and.l_ri==.false.) then 
!!!!!     write(6,*)'single precision integrals'
         call run_rhf_s(icore, isympopa, isympopb, isymvrta, isymvrtb, isymmo, &
                                      nvmax, nomax, nmo) 
         elseif(prec==4.and.l_ri==.true.) then 
         call run_rhf_ri_s(icore, isympopa, isympopb, isymvrta, isymvrtb, isymmo, &
                                      nvmax, nomax, nmo) 
         endif
     elseif(iop_ref==1) then 
!         write(6,*)'isympopb', isympopb(1,2)
!         write(6,*)'isymvrtb', isymvrtb
         if(prec==8.and.l_ri==.False.) then 
            call run_uhf(icore, isympopa, isympopb, isymvrta,isymvrtb, isymmo_u, &
                                      nvmax, nomax, nmo)
         elseif(prec==4.and.l_ri==.False.) then 
            call run_uhf_s(icore, isympopa, isympopb, isymvrta,isymvrtb, isymmo_u, &
                                      nvmax, nomax, nmo)
         elseif(prec==8.and.l_ri==.True.) then 
            call run_uhf_ri(icore, isympopa, isympopb, isymvrta, isymvrtb, isymmo_u, &
                                      nvmax, nomax, nmo) 
         elseif(prec==4.and.l_ri==.True.) then 
            call run_uhf_ri_s(icore, isympopa, isympopb, isymvrta, isymvrtb, isymmo_u, &
                                      nvmax, nomax, nmo) 
         else
            write(6,*)'error in CCSOINP, iop_ref, prec, l_ri'
         endif
     endif

     call flush(6)
     call walltimeso(rt2)
     write(6,*)'Time for cc_pyscf----', rt2-rt1

      !=========for old ccsdso
     i0 = 1
     if(prec==8.and.eri_cfour==1) then
     call clmoioso0(1, icore, pop(1,1), pop(1,2), 0, &
                             pop(1,1), vrt(1,2), 1, 'IJKAAA__')
     call clmoioso0(1, icore, pop(1,1), pop(1,2), 1,&
                             vrt(1,1), pop(1,2), 1, 'IJAKAB__')
     call getallso(icore, 1, 'IJKAAB__')
     i1 = i0 + idsymsoc(1, popa, popa, 1, popa, vrta, 1)*isd
     call iassymso(icore, icore(i1), 1, popa, popa, popa, vrta, 5)
     call putallso(icore(i1), 1, 'IJKAAA__')
!   <j|Ak> = <jI|kA> 
     call isymtrso(icore, icore(i1), popa, popa, 1, 1,&
                               popa, vrta, 1, 0, 1)
     call isymtrso(icore(i1), icore, popa, popa, 1, 0,&
                               popa, vrta, 1, 1, 1)
     call putallso(icore, 1, 'IJAKAB__')
!     write(6,*)'finished ijka'
     if(eri_level==1.or.eri_level==0) then 
       if(.not.l_ri) then
       call clmoioso0(1, icore, vrt(1,1), vrt(1,2), 0,&
                          vrt(1,1), pop(1,2), 1, 'ABCIAA__')
       endif
       if(labic) then
       call clmoioso0(1, icore, vrt(1,1), vrt(1,2), 1,&
                          pop(1,1), vrt(1,2), 1, 'ABICAB__')
       endif
     ioffs8(1) = 1
       do irrepci = 1, nirrep-1
          irrepab = irrepci 
          numci = irpdso(irrepci, vrt, pop, 1)
          numab = irpdso(irrepab, vrt, vrt, 1)
          ioffs8(irrepci+1) = ioffs8(irrepci)+numab*numci
       enddo
     i0 = 1
     i2 = i0 + isd*idsymsoc(1, vrta, vrta, 1, vrta, popa, 1)*2
     write(6,*)'get abciab in small core'
     if((.not.l_ri).or.labic) then
        do irrepi = 1, nirrep
           numi = popa(irrepi)
           nsize = idsymsoc(irrepi, zorb, vrta, 1, vrta, vrta, 1)
           i1 = i0 + isd*idsymsoc(irrepi, zorb, vrta, 1, vrta, vrta, 1) 
           i2 = i1 + isd*idsymsoc(irrepi, zorb, vrta, 1, vrta, vrta, 1) 
           do i = 1, popa(irrepi)
              call getpqrso(icore(i0), i, irrepi, 1, 'ABCIAB__')
              if (labic) then
              call isymtrso(icore(i0), icore(i2), vrt, vrt, 1, 1, vrt,&
                                    zorb, 1, 0, irrepi)
!              write(6,*)'norm of abc', i, nrm2so(nsize, icore, 1)
              ioffabc(1) = 1 
              do irrepc = 1, nirrep
                 irrepci = dirprd(irrepc, irrepi) 
                 irrepab = dirprd(irrepci, 1)
                 numab = irpdso(irrepab, vrta, vrta, 1) 
                 numc = vrta(irrepc)
                 ioffabc(irrepc+1) = ioffabc(irrepc) + numab*numc*isd 
                 ioffr0 = ioffabc(irrepc)
                 ioffs0 = (isymoffso(irrepc, irrepci, pop, vrt,1)-1)+1
!                 write(6,*)'ioffs', ioffabc
                 do c = 1, numc
                    ioffr = ioffr0 + (c-1)*numab*isd
                    ioffs = ioffs0 + ((c-1)*numi+i-1)
!                   ioffs2 = (ioffs0-1 + (c-1)*numi+i-1)*numab &
!                             + ioffs8(irrepci)
!                   write(6,*)'ioffs', ioffr, ioffs 
                    call putlstso(icore(i2-1+ioffr), ioffs, 1,&
                            irrepci, 1, 'ABICAB__')
!                   ioffs2 = ((c-1)*numi + i-1)*numab + 1 
!                   call myicopyso(icore(ioffr),
!                           icore(i2-1+ioffs2), numab)
                 enddo
               enddo
!              write(6,*)'norm of abc', i, nrm2so(nsize, icore, 1)
              endif
              call iassymso(icore(i0), icore(i1), irrepi, vrta, vrta,&
                             vrta, zorb, 5)
              call putpqrso(icore(i1), i, irrepi, 1, 'ABCIAA__')
           enddo
        enddo
     endif
     nsize = idsymsoc(1,vrta,vrtb,1 ,vrta,popb, 1)
!    call getallso(icore, 1, 'ABICAB__')
!    write(6,*)'norm of ABICAB  ',nsize, nrm2so(nsize, icore, 1)
!    write(6,*)'norm of ABICAB  ',nsize, nrm2so(nsize, icore(i2), 1)
!    call getallso(icore, 1, 'ABCIAB__')
!    write(6,*)'norm of ABCIAB  ',nsize, nrm2so(nsize, icore, 1)
!    call getallso(icore, 1, 'ABCIAA__')
!    nsize = idsymsoc(1,vrta,vrtb,0 ,vrta,popb, 1)
!    write(6,*)'norm of ABCIAA  ',nsize, nrm2so(nsize, icore, 1)

!      call getallso(icore, 1, 'ABCIAB__')
!      i1 = i0 + idsymsoc(1, vrta, vrta, 1, vrta, popa, 1)*isd
!      call iassymso(icore, icore(i1), 1, vrta, vrta, vrta, popa, 5)
!      call putallso(icore(i1), 1, 'ABCIAA__')

!      call isymtrso(icore, icore(i1), vrta, vrta, 1, 1,
!                              vrta, popa, 1, 0, 1)
!      call isymtrso(icore(i1), icore, vrta, vrta, 1, 0,
!                              vrta, popa, 1, 1, 1)
!      call putallso(icore, 1, 'ABICAB__')
     endif
!    call getallso(icore, 1, 'ABICAB__')
!    nsize = idsymsoc(1,vrta,vrtb,1 ,vrta,popb, 1)
!    write(6,*)'norm of ABICAB  ',nsize, nrm2so(nsize, icore, 1)
!    call getallso(icore, 1, 'ABCIAA__')
!    nsize = idsymsoc(1,vrta,vrtb,0 ,vrta,popb, 1)
!    write(6,*)'norm of ABCIAA  ',nsize, nrm2so(nsize, icore, 1)

!    write(6,*)'finished abci'

     call clmoioso0(1, icore, vrt(1,1), pop(1,2), 1,&
                           vrt(1,1), pop(1,2), 1, 'ABIJAB2_')
     call clmoioso0(1, icore, vrt(1,1), vrt(1,2), 0,&
                           pop(1,1), pop(1,2), 0, 'ABIJAA__')
!     write(6,*)'finished abij'
     call getallso(icore, 1, 'ABIJAB__')
     i1 = i0 + idsymsoc(1, vrta, vrta, 1, pop, popa, 1)*isd
     call sstgenso(icore, icore(i1), nsize, vrta, vrta, popa, popa,&
                   icore, 1, '1423')
!     write(6,*)'finished abij'
     call putallso(icore(i1), 1, 'ABIJAB2_')

!     write(6,*)'finished abij'
     call iassymso(icore, icore(i1), 1, vrta, vrta, popa, popa, 5)
     call icompso(icore(i1), icore, vrta, vrta,0,1, popa, popa,1,0,1)
     call putallso(icore, 1, 'ABIJAA__')
!     write(6,*)'finished abij'
!<AB|IJ> to <AI,BJ>
     call clmoioso0(1, icore, vrt(1,1), pop(1,2), 1,&
                           vrt(1,1), pop(1,2), 1, 'ABIJAA2_')
     call iexpso(icore(i1), icore, vrta, vrta,0,1, popa, popa,1,0,1)
     call sstgenso(icore, icore(i1), nsize, vrta, vrta, popa, popa,&
                   icore, 1, '1324')
     call putallso(icore(i1), 1, 'ABIJAA2_')
!     write(6,*)'finished abij'

     call clmoioso0(1, icore, pop(1,1), pop(1,2), 0,&
                           pop(1,1), pop(1,2), 0, 'IJKLAA__')
     call getallso(icore, 1, 'IJKLAB__')
     i1 = i0 + idsymsoc(1, popa, popa, 1, pop, popa, 1)*isd
     call iassymso(icore, icore(i1), 1, popa, popa, popa, popa, 5)
     call icompso(icore(i1), icore, popa, popa,0,1, popa, popa,1,0,1)
     call putallso(icore, 1, 'IJKLAA__')
!     write(6,*)'finished ijkl'

     call getallso(icore, 1, 'BIAJABAB')
     nsize = idsymsoc(1,vrta,popb,1, vrta, popb,1)
     call scalso(nsize, -1.d0, icore, 1)
     call putallso(icore, 1, 'BIAJABAB')

     call getallso(icore, 1, 'BIAJAAAA')
     nsize = idsymsoc(1,vrtb,popa, 1,vrtb, popa,1)
     call scalso(nsize, -1.d0, icore, 1)
     call putallso(icore, 1, 'BIAJAAAA')

!
     if(l_ri==.true..and.eri_level==0) then
     call clmoioso0(1, icore, vrt, vrt, 0,&
                           vrt, vrt, 0, 'ABCDAA__')
     i0 = 1 
        do irrepd = 1, nirrep 
           i1 = i0 + isd*idsymsoc(irrepd, zorb, vrta, 1, vrta, vrta, 1) 
           do d = 1, vrta(irrepd) 
              call getpqrso(icore(i0), d, irrepd, 1, 'ABCDAB__')
              call iassymso(icore(i0), icore(i1), irrepd, vrta, vrta, &
                             vrta, zorb, 5)
              call putpqrso(icore(i1), d, irrepd, 1, 'ABCDAA__') 
           enddo
        enddo
     endif
     call getallso(icore, 1, 'IJKLAA__')
     nsize = idsymsoc(1, popa, popb, 0, popa, popb, 0)
     write(6,*)'norm of IJKLAA  ',nsize, nrm2so(nsize, icore, 1)

     call getallso(icore, 1, 'IJKLAB__')
     nsize = idsymsoc(1, popa, popb, 1, popa, popb, 1)
     write(6,*)'norm of IJKLAB  ',nsize, nrm2so(nsize, icore, 1)

     call getallso(icore, 1, 'IJKAAA__')
     nsize = idsymsoc(1, popa, popb, 0, popa, vrtb, 1)
     write(6,*)'norm of IJKAAA  ',nsize, nrm2so(nsize, icore, 1)

     call getallso(icore, 1, 'IJKAAB__')
     nsize = idsymsoc(1, popa, popb, 1, popa, vrtb, 1)
     write(6,*)'norm of IJKAAB  ',nsize, nrm2so(nsize, icore, 1)

     call getallso(icore, 1, 'IJAKAB__')
     nsize = idsymsoc(1, popa, popb, 1, popa, vrtb, 1)
     write(6,*)'norm of IJAKAB  ',nsize, nrm2so(nsize, icore, 1)

     call getallso(icore, 1, 'ABIJAA__')
     nsize = idsymsoc(1, vrta,vrtb,0,popa,popb, 0)
     write(6,*)'norm of ABIJAA  ',nsize, nrm2so(nsize, icore, 1)

     call getallso(icore, 1, 'ABIJAB__')
     nsize = idsymsoc(1, vrta,vrtb,1,popa,popb, 1)
     write(6,*)'norm of ABIJAB  ',nsize, nrm2so(nsize, icore, 1)

     call getallso(icore, 1, 'BIAJAABB')
     nsize = idsymsoc(1, vrta,popa,1,vrtb,popb, 1)
     write(6,*)'norm of BIAJAABB',nsize, nrm2so(nsize, icore, 1)

     call getallso(icore, 1, 'BIAJABAB')
     nsize = idsymsoc(1,vrta,popb,1, vrta, popb,1)
     write(6,*)'norm of BIAJABAB',nsize, nrm2so(nsize, icore, 1)

     call getallso(icore, 1, 'BIAJAAAA')
     nsize = idsymsoc(1,vrtb,popa, 1,vrtb, popa,1)
     write(6,*)'norm of BIAJAAAA',nsize, nrm2so(nsize, icore, 1)

!    call getallso(icore, 1, 'ABCIAA__')
!    nsize = idsymsoc(1,vrta,vrtb,0 ,vrta,popb, 1)
!    write(6,*)'norm of ABCIAA  ',nsize, nrm2so(nsize, icore, 1)

!    call getallso(icore, 1, 'ABCIAB__')
!    nsize = idsymsoc(1,vrta,vrtb,1 ,vrta,popb, 1)
!    write(6,*)'norm of ABCIAB  ',nsize, nrm2so(nsize, icore, 1)

!    call getallso(icore, 1, 'ABICAB__')
!    nsize = idsymsoc(1,vrta,vrtb,1 ,vrta,popb, 1)
!    write(6,*)'norm of ABICAB  ',nsize, nrm2so(nsize, icore, 1)

!    if(ldebug.and.eri_level==0) then
!    call getallso(icore, 1, 'ABCDAA__')
!    nsize = idsymsoc(1,vrta,vrtb,0 ,vrta,vrtb, 0)
!    write(6,*)'norm of ABCDAA  ', nrm2so(nsize, icore, 1)

!    call getallso(icore, 1, 'ABCDAB__')
!    nsize = idsymsoc(1,vrta,vrtb,1 ,vrta,vrtb, 1)
!    write(6,*)'norm of ABCDAB  ', nrm2so(nsize, icore, 1)

!    call getallso(icore, 1, 'ABCDAB2_')
!    nsize = idsymsoc(1,vrta,vrtb,2 ,vrta,vrtb, 2)
!    write(6,*)'norm of ABCDAB2 ', nrm2so(nsize, icore, 1)
!    endif
          
     i0 = 1
     i1 = i0 + nva**2*isd
     i2 = i1 + nva**2*isd
     i3 = i2 + nva**2*isd
     i4 = i3 + nva**2*isd
     i5 = i4 + nva**2*isd
     i6 = i5 + nva**2*isd
     i7 = i6 + nva**2*isd
     i8 = i7 + nva**2*isd
     i9 = i8 + nva**2*isd
     i10= i9 + nva**2*isd
     i11= i10 + nva**2*isd
     i12= i11 + nva**2*isd
     call SAVEINTEGRALS(NOA, NVA, ESCF, EA_resort, icore(i12), icrsiz,&
                        icore(i0), icore(i1),icore(i2), icore(i3),&
                        icore(i4), icore(i5),&
                        icore(i6), icore(i7), icore(i8),&
                        icore(i9), icore(i10), icore(i11))
      endif 
      write(6,*) 'finish cc_pyscf'
 
      END

